Electronic stopping power in insulators from first principles 
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Using time-dependent density-functional theory we calculate from first principles the rate of energy 
transfer from a moving proton or antiproton to the electrons of an insulating material, LiF. The 
behavior of the electronic stopping power versus projectile velocity displays an effective threshold 
velocity of ~ 0.2 a.u. for the proton, consistent with recent experimental observations, and also for 
the antiproton. The calculated proton/antiproton stopping-power ratio is 2.4 at velocities slightly 
above the threshold {v ~ 0.4 a.u.), as compared to the experimental value of 2.1. The projectile 
energy loss mechanism is observed to be stationary and extremely local. 

PACS numbers: 61.85.+p; 61.80.Jh; 61. 82. Ms; 34.50. Bw; 71.15.Pd 
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The interaction of charged particles with matter has 
been a subject of extensive research since the discovery 
of subatomic particles Ions moving through solids 
gradually transmit their kinetic energy to electronic ex- 
citations of the host and deposit it along their path. 
The maximum of this deposited energy is the so-called 
Bragg peak and it occurs shortly before the particle stops. 
Hence the importance of studying the electronic energy 
loss of slow ions (-with velocities belo'w the Bohr velocity) 
travelling through solids. For metals, the understand- 
ing of this problem has been steadily progressing over 
the years [Ij. For insulators, however, experimental re- 
sults remain unexplained even for simple systems. This 
is particularly true at low velocities (the threshold effect) 
@) Si I where the contribution from nuclear collisions con- 
ceals the electronic stopping [Jj. 

A detailed quantitative knowledge of these processes 
is required to understand the damage produced in mate- 
rials when exposed to radiation. For ceramic materials 
devised for the encapsulation of nuclear waste 0] the pre- 
diction of durability over extremely long times is crucial. 
Radiation-damage simulations performed to date 0, 01 
rely on empirical force fields obtained from fits to low- 
energy properties. The actual interatomic forces could be 
enormously altered, however, by the local electron heat- 
ing produced by the electronic stopping. 

In the semiclassical formalism, the electronic energy 
loss rate is given by the response of the system to the 
external potential, ^ = — Zv • E'"*^, where Z and v are 
the charge and velocity of the projectile, and E™"* is the 
induced electric field in the target material. Hence, to 
first order, linear response theory can be used to give 
the stopping power (SP) -the energy lost per unit path 
length- in terms of the dynamical dielectric function of 
the material. For metals this formalism shows that the 
SP is hnear with small velocities, ^ ^ {Zefv [1, 



reflecting that no minimum energy is required to excite 
electron-hole pairs. It must be remembered, however, 
that the dielectric approach is not valid at low ion veloc- 
ities, and non-linear effects cannot be neglected 

A different behavior is expected in wide band-gap insu- 
lators, given the finite energy required to excite electrons. 
For protons moving through noble gases , which also 
display a large minimum excitation energy, a velocity 



threshold has been observed experimentally 11| , and ex- 
plained in terms of the quantization of energy levels and 
charge exchanges. However, for solid insulators the situ- 
ation is unclear. No threshold effects were originally ob- 
served in AI2O3, Si02, or LiF, and the hnear dependence 
dE/dx oc V was observed down to velocities of about 0.3 

a.u. iEm. 

For protons under grazing incidence in LiF(OOl), and 
below ~ 0.2 a.u. a threshold behavior was reported [l4| . 
Under these conditions the proton does not penetrate 
the solid, and charge exchange is identified as the domi- 
nant mechanism for electronic stopping, with local elec- 
tron capture from F~ ions, giving rise to H*^ and H^. 
More recent experiments on thin LiF films show an ap- 
parent velocity threshold near 0.1 a.u. The differ- 
ent experimental setup (transmission) suggests a differ- 
ent physical stopping mechanism, based on electron-hole 
pair excitations. A threshold behavior is expected, and 
has been qualitatively predicted from linear-response and 
from perturbation theory calculations. These approxima- 
tions, however, grossly underestimate the SP at low ve- 
locities, thus exaggerating the threshold [3]. A theoret- 
ical description beyond these approximations is needed. 

In this letter we present a first-principles approach 
to the non-perturbative study of realistic solid-ion in- 
teractions. We follow the explicit time evolution of the 
electronic states of the host crystal as an external par- 
ticle propagates through the system, by means of time- 



2 




V, = 0.387 



O 

o 



V =0.096 



20 



15 



10 



10 12 14 



x(A) 



FIG. 1: Total electronic energy as a function of displacement 
for several proton velocities. 



dependent density- functional theory (TD-DFT) [l^] ■ We 
use here this scheme to understand the SP threshold ef- 
fect and stopping mechanisms in LiF, a well-studied char- 
acteristic insulator, finding reasonable agreements with 
measured data. Most importantly, however, this study 
sets the scene for a promising line of theoretical simula- 
tions, assesses its possibilities and offers new insights into 
the stopping process. 

The energy transmitted to the electrons from a 
constant-velocity moving ion is monitored. Total energy 
is thus not conserved since the energy loss of the projec- 
tile will be neglected (its large mass ensures a negligibly 
small decline in its velocity on the time scales of the simu- 
lations). As the center of the charged particle, Tc, moves 
with constant velocity, = ro+v-t, the time-dependent 
Kohn-Sham (KS) equation [16] defines the dynamics of 
effective single-particle states (and thereby the electronic 
density and energy) under the external potential gener- 
ated by the projectile and the crystal of Li and F nuclei. 

The calculations were done using the Siesta ab initio 
method [itI . [l^ , in its time-evolving TD-DFT implemen- 
tation 'l9'|, using the instantaneous local density approx- 
imation (LDA) to exchange and correlation [20|. The Is 
core electron pair of F was replaced by a norm-conserving 
pseudopotential [ilj in the fully nonlocal form. The Is 
electrons of Li were explicitly included in the calculation 
as pseudized valence electrons. A double- ^ polarized ba- 
sis was considered for the valence electrons (sirigle-C for 
the Is of Li), with an energy shift of 100 meV [18|. The 
grid cutoff [18[ for integration was 118 Ry. The lattice pa- 
rameter obtained for bulk LiF is 3.98 A, slightly smaller 
than the experimental value 4.03 A, as expected for LDA. 
The projectile (proton, p, or antiproton, p) was described 
by the bare ±i potential It was not pseudized in 

order to treat p and p on the same footing. 

Periodic boundary conditions were used throughout. 
The supercell size was chosen so as to minimize the 
spurious effects of the repetition while keeping manage- 



able computational demands. After convergence tests, 
a 4x4x4 supercell of 128 atoms was selected, such that 
a particle moves ^^llA down a [110] channel before re- 
entering the box. The F-point approximation was used 
for integrations in k. The projectile was initially put in 
the centre of a crystal cage and the time-independent 
DFT solution was obtained to define the initial state for 
the subsequent evolution. It was then moved with con- 
stant velocity along a [110] channel. The time-dependent 
KS equation was then solved numerically by discretizing 
time and applying the Crank-Nicholson algorithm as de- 
scribed in ref. 19]. Using a time step of 1 attosecond, 
the wavefunctions were then propagated for several fem- 
toseconds. The electronic total energy was recorded as a 
function of time and the SP, ^ , was extracted. 

Fig. [T] shows the electronic total energy of the system 
as a function of displacement for different velocities of 
the projectile. At low velocities the adiabatic behavior 
is recovered, with no net energy transfer, just the ex- 
pected oscillation of the total energy with the position 
of the projectile in the crystal. At higher velocities the 
oscillations are superimposed to an underlying energy in- 
crease with time. After a remarkably short transient pe- 
riod (around 0.3 fs), the energy increase stabilizes to an 
essentially stationary regime, in which the energy differ- 
ence between consecutive equivalent lattice positions of 
the projectile remains constant. The SP is then extracted 
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FIG. 2: Calculated electronic stopping power -^j;^ as a function 
of particle velocity v for p (empty circles) and p (full circles). 
Crosses indicate calculations for p including additional basis 
orbitals along the projectile's path (see text). The inset shows 
the excess energy acquired by a octahedral LigF^"'" cluster on 
the passage of a p on an analogous trajectory to the solid's. 
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from the average slope of this stationary regime. Stop- 
ping power results for p and p are presented in Fig. [2] A 
threshold effect at low velocities is apparent in the fig- 
ure, in agreement with recent experiments (sj, unlike the 
linear behavior observed for metals [l| , but much smaller 
than that predicted for insulators by linear response the- 
ory [l^. The values obtained in this work (around 0.2 
a.u.) are consistent with experiments 0, The differ- 
ence between p and p is also apparent, in contrast with 
the invariance under charge-sign change expected from 
perturbative treatments. Experiments show that the en- 
ergy loss for protons is approximately twice as high as 
for antiprotons [3] , consequence of the different kind of 
screening (one attracts electrons while the other attracts 
holes). In the neighborhood oi v — 0.4 a.u. (slightly 
above the threshold) we obtain a SP ratio of ~ 2.4, that 
compares well with the experimental value of 2.1 [13i] . 

Since a charge dependence is expected from linear 
response, the SP over q"^ is plotted in Fig. [3] for both 
projectiles. The smooth behavior at the origin indicates 
that there are no substantial biases in the way the stop- 
ping mechanism is described for p and p. The slope near 
the origin corresponds to the dependence, while the 
bending corresponds to q^ and higher terms. The dis- 
played behavior, and particularly the sign of both terms, 
is similar to that observed for metals [22i |. 

There is a discrepancy between measured and com- 
puted values for the SP for velocities above the threshold. 
The experiments show four times higher SP [3| . A factor 
of two is accounted for by the known relation between 
channelling conditions and the average over random tra- 



jectories [23|, |2J]. Technical reasons account for a part 
of the remaining discrepancy. The inclusion of basis or- 
bitals along the projectile's path (a sp single-C set every 
0.5 A) increases the SP from 1.6 eV/A to 2.8 eV/A for 
V — 0.46 a.u. [3^ as shown in Fig. [21 At a more funda- 
mental level, discrepancies are expected from the errors 
in electronic spectra around the band gap given by in- 
stantaneous LDA in TD-DFT (even if it includes a RPA 
correction on KS eigenvalues) [25|. 
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FIG. 3: Dependence of with the charge q of the external 



We now analyze the locality of the energy absorption in 
terms of MuUiken charges [26| . Although they are known 
to be unreliable in their explicit dependence on the basis 
set, their changes for fixed basis set are much more in- 
dicative. By comparing the atomic charges obtained dy- 
namically and adiabatically (a static self-consistent cal- 
culation with the external potential at the same instan- 
taneous position) , insights are obtained into the stopping 
process. Fig.|4]shows this comparison for the proton case, 
where the following is observed: (z) The screening of the 
proton is enhanced in the dynamic case (our basis set 
describes this screening by a transfer of charge from F to 
Li), due to the increased polarizability of the host at the 
frequencies explored by the moving ion. (ii) A delay is 
apparent in the dynamical screening, which is the main 
responsible of the dissipation, (in) The electronic exci- 
tation process is extremely local: the changes observed 
in atoms closest to the trajectory are much larger than 
any other, and the dynamical screening effect is only no- 
ticeable when the projectile is very close to the ion. 

The locality of the energy transfer is confirmed when 
calculating the energy absorbed by a small cluster of LiF 
(LigF^+). The inset of Fig. [2] shows a striking similarity 
in the overall w-dependence of the energy absorbed by the 
cluster and the SP in the solid. If we took the effective 
path length in the cluster as 1.4 A (the number of va- 
lence electrons in the cluster corresponds to one formula 
unit), both SPs would be indeed very similar 3J|. In 
the solid, the energy accumulated along the path would 
then diffuse away at longer time scales (corresponding to 
the effective band width), defining a wake. A tentative 
definition of local energy [3^ shows well differentiated 
time scales for the excitation by the projectile, on one 
hand, and the ensuing out-diffusion, on the other (not 
shown). The short-ranged initial excitation can be ratio- 
nalized in terms of the electronic localization length scale 
relevant to dielectric response 27|, which is expected to 
be very short for LiF. Similar locality in the energy loss 
mechanism has been recently found in a confined two- 
dimensional electron gas, however [28| . 

A characterization of the charge state of the projec- 
tile 



projectile, for v = 0.46 a.u. The stopping power is in eV/A, 
and q in electron charge units. 



30| has not been attempted. The present calcu- 
lations allow for the establishment of charge states of any 
kind within the constraints that the electronic charge and 
momentum are conserved. It is clear that the mid-gap 
state travels with the proton, and that it becomes partly 
populated. It is not clear, however, that the charge as- 
sociated to that state represents a meaningful definition 
of the charge state of the projectile, since whether the 
screening charge builds up as it passes or travels with it 
is not determined by the population of that state. This 
effect will be explored in further works. 

In conclusion, we have presented a general approach 
for the non-perturbative first-principles study of the elec- 
tronic stopping power in solids. New insights into the 
electronic SP in insulators have been provided for pro- 
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FIG. 4: Electron population on Li atoms close to the proton 
trajectory versus position of the projectile along it. Adiabatic 
(dynamic) charges are shown by circles (dashed lines). The 
position of the maximum in the adiabatic curve occurs when 
the projectile is closest to the atom (~1.2A). The almost con- 
stant line corresponds to a second closest Li atom, 2.3A away 
from the trajectory. 



tons and antiprotons in LiF. 
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